Section 1: Logistic regression

We will analyse the data collected by Jones (Unpublished BSc dissertation, University of Southampton, 1975). The aim of the study was to define if the probability of having Bronchitis is influenced by smoking and/or pollution.

The data are stored under data/Bronchitis.csv and contains information on 212 participants.

Section 1.1: importation and descriptive analysis

Lets starts by

importing the data set Bronchitis with the function read.csv()
displaying bron (a dichotomous variable which equals 1 for participants having bronchitis and 0 otherwise) as a function of cigs, the number of cigarettes smoked daily.

Bronchitis = read.csv("data/Bronchitis.csv",header=TRUE)
plot(Bronchitis$cigs,Bronchitis$bron,col="blue4",
        ylab = "Absence/Presense of Bronchitis", xlab = "Daily number of cigarettes")
abline(h=c(0,1),col="light blue")

Section 1.2: Model fit

Lets

fit a logistic model by means the function glm() and by means of the function gamlss() of the library gamlss.
display and analyse the results of the glm function : Use the function summary() to display the results of an R object of class glm.

fit.glm = glm(bron~cigs,data=Bronchitis,family=binomial)

library(gamlss)
fit.gamlss = gamlss(bron~cigs,data=Bronchitis,family=BI)

## GAMLSS-RS iteration 1: Global Deviance = 181.7072 
## GAMLSS-RS iteration 2: Global Deviance = 181.7072

summary(fit.glm)

## 
## Call:
## glm(formula = bron ~ cigs, family = binomial, data = Bronchitis)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4418  -0.5472  -0.4653  -0.4405   2.1822  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.2840     0.2731  -8.365  < 2e-16 ***
## cigs          0.2094     0.0376   5.567 2.59e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 221.78  on 211  degrees of freedom
## Residual deviance: 181.71  on 210  degrees of freedom
## AIC: 185.71
## 
## Number of Fisher Scoring iterations: 4

Let’s now define the estimated probability of having bronchitis for any number of daily smoked cigarette and display the corresponding logistic curve on a plot:

plot(Bronchitis$cigs,Bronchitis$bron,col="blue4",
        ylab = "Absence/Presense of Bronchitis", xlab = "Daily number of cigarettes")
abline(h=c(0,1),col="light blue")

axe.x = seq(0,40,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])/(1+exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2]))
lines(axe.x,f.x,col="pink2",lwd=2)

Section 1.3: Model selection

As for linear models, model selection may be done by means of the function anova() used on the glm object of interest.

anova(fit.glm,test="LRT")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: bron
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev Pr(>Chi)    
## NULL                   211     221.78             
## cigs  1    40.07       210     181.71 2.45e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Section 1.3: Model check

Lets assess is the model fit seems satisfactory by means

of the analysis of deviance residuals (function plot() on an object of class glm,
of the analysis of randomised normalised quantile residuals (function plot() on an object of class gamlss,

# deviance
par(mfrow=c(2,2),mar=c(3,5,3,0))
plot(fit.glm)

# randomised normalised quantile residuals
plot(gamlss(bron~cigs,data=Bronchitis,family=BI))

## GAMLSS-RS iteration 1: Global Deviance = 181.7072 
## GAMLSS-RS iteration 2: Global Deviance = 181.7072

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.05227133 
##                        variance   =  0.9386897 
##                coef. of skewness  =  0.02624943 
##                coef. of kurtosis  =  2.950084 
## Filliben correlation coefficient  =  0.9984055 
## ******************************************************************

Section 1.4: Fun

# long format:
long = data.frame(mi = rep(c("MI","No MI"),c(104+189,11037+11034)),
                  treatment = rep(c("Aspirin","Placebo","Aspirin","Placebo"),c(104,189,11037,11034)))
# short format: 2 by 2 table
table2by2 = table(long$treatment,long$mi)

#
chisq.test(table2by2)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table2by2
## X-squared = 23.782, df = 1, p-value = 1.079e-06

prop.test(table2by2[,"MI"],apply(table2by2,1,sum))

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  table2by2[, "MI"] out of apply(table2by2, 1, sum)
## X-squared = 23.782, df = 1, p-value = 1.079e-06
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.010570828 -0.004440228
## sample estimates:
##      prop 1      prop 2 
## 0.009334889 0.016840417

summary(glm(I(mi=="MI")~treatment,data=long,family="binomial"))

## 
## Call:
## glm(formula = I(mi == "MI") ~ treatment, family = "binomial", 
##     data = long)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.1843  -0.1843  -0.1370  -0.1370   3.0574  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -4.66462    0.09852 -47.348  < 2e-16 ***
## treatmentPlacebo  0.59763    0.12283   4.865 1.14e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3122.5  on 22363  degrees of freedom
## Residual deviance: 3097.8  on 22362  degrees of freedom
## AIC: 3101.8
## 
## Number of Fisher Scoring iterations: 7

Section 2: Poisson regression

The dataset students.csv shows the number of high school students diagnosed with an infectious disease for each day from the initial disease outbreak.

Section 2.2: Importation

Lets

import the dataset by means of the function read.csv()
display the daily number of students diagnosed with the disease (variable cases) as a function of the days since the outbreak (variable day).

students = read.csv("data/students.csv",header=TRUE)
plot(students$day,students$cases,col="blue4",
        ylab = "Number of diagnosed students", xlab = "Days since initial outbreak")
abline(h=c(0),col="light blue")

Section 2.2: Model fit

Lets

fit a poisson model by means the function glm() and by means of the function gamlss() of the library gamlss.
display and analyse the results of the glm function : Use the function summary() to display the results of an R object of class glm.

fit.glm = glm(cases~day,data=students,family=poisson)

library(gamlss)
fit.gamlss = gamlss(cases~day,data=students,,family=PO)

## GAMLSS-RS iteration 1: Global Deviance = 389.1082 
## GAMLSS-RS iteration 2: Global Deviance = 389.1082

summary(fit.glm)

## 
## Call:
## glm(formula = cases ~ day, family = poisson, data = students)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.00482  -0.85719  -0.09331   0.63969   1.73696  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.990235   0.083935   23.71   <2e-16 ***
## day         -0.017463   0.001727  -10.11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 215.36  on 108  degrees of freedom
## Residual deviance: 101.17  on 107  degrees of freedom
## AIC: 393.11
## 
## Number of Fisher Scoring iterations: 5

Let’s now define the estimated probability of having bronchitis for any number of daily smoked cigarette and display the corresponding logistic curve on a plot:

plot(students$day,students$cases,col="blue4",
        ylab = "Number of diagnosed students", xlab = "Days since initial outbreak")
abline(h=c(0),col="light blue")

axe.x = seq(0,120,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])
lines(axe.x,f.x,col="pink2",lwd=2)

Section 2.3: Model selection

As for linear models, model selection may be done by means of the function anova() used on the glm object of interest.

anova(fit.glm,test="LRT")

## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: cases
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                   108     215.36              
## day   1   114.18       107     101.17 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Section 2.3: Model check

Lets assess is the model fit seems satisfactory by means

of the analysis of deviance residuals (function plot() on an object of class glm,
of the analysis of randomised normalised quantile residuals (function plot() on an object of class gamlss,

# deviance
par(mfrow=c(2,2),mar=c(3,5,3,0))
plot(fit.glm)

# randomised normalised quantile residuals
plot(fit.gamlss)

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  0.01029368 
##                        variance   =  0.895386 
##                coef. of skewness  =  -0.444238 
##                coef. of kurtosis  =  2.638977 
## Filliben correlation coefficient  =  0.9881143 
## ******************************************************************

Section 6: Practicals

(i) Bronchitis.csv

Analyse further the Bronchitis data of Jones (1975) by

first investigating if the probability of having bronchitis also depends on pollution (variable poll),
second investigating if there is an interaction between the variables cigs and poll.

Lets plot the data first.

Bronchitis = read.csv("data/Bronchitis.csv",header=TRUE)
# plot
plot(Bronchitis$poll,Bronchitis$bron,col="blue4",
        ylab = "Absence/Presense of Bronchitis", xlab = "Pollution level")
abline(h=c(0,1),col="light blue")

No obvious relationship between pollution and bronchitis is visible by means of this plot.

Lets fit a model assuming that the probability of getting bronchitis is a function of the pollution level.

# fit1:
fit.glm = glm(bron~poll,data=Bronchitis,family=binomial)
summary(fit.glm)

## 
## Call:
## glm(formula = bron ~ poll, family = binomial, data = Bronchitis)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1018  -0.7254  -0.6004  -0.4944   2.0796  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -8.65578    2.59070  -3.341 0.000834 ***
## poll         0.12482    0.04341   2.876 0.004032 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 221.78  on 211  degrees of freedom
## Residual deviance: 213.26  on 210  degrees of freedom
## AIC: 217.26
## 
## Number of Fisher Scoring iterations: 4

The intercept of the previous fit allows to define the probability of getting bronchitis when the level of pollution equals 0. As a zero level of pollution is (i) out of range (ii) not a realistic value, we will

create the variable poll_centered defined as the pollution level minus the mean so that the intercept corresponds to the probability of getting bronchitis for an average pollution level in Cardiff,
refit the model

# fit2:
Bronchitis$poll_centered = Bronchitis$poll-mean(Bronchitis$poll)
fit.glm = glm(bron~poll_centered,data=Bronchitis,family=binomial)
library(gamlss)
fit.gamlss = gamlss(bron~poll_centered,data=Bronchitis,family=BI)

## GAMLSS-RS iteration 1: Global Deviance = 213.2586 
## GAMLSS-RS iteration 2: Global Deviance = 213.2586

summary(fit.glm)

## 
## Call:
## glm(formula = bron ~ poll_centered, family = binomial, data = Bronchitis)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.1018  -0.7254  -0.6004  -0.4944   2.0796  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.34808    0.17576  -7.670 1.72e-14 ***
## poll_centered  0.12482    0.04341   2.876  0.00403 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 221.78  on 211  degrees of freedom
## Residual deviance: 213.26  on 210  degrees of freedom
## AIC: 217.26
## 
## Number of Fisher Scoring iterations: 4

Lets

perform a model check but plotting the randomised quantile residuals of gamlss a few times
plot the fitted probabilities

# model check
plot(gamlss(bron~poll_centered,data=Bronchitis,family=BI))

## GAMLSS-RS iteration 1: Global Deviance = 213.2586 
## GAMLSS-RS iteration 2: Global Deviance = 213.2586

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  0.02469888 
##                        variance   =  0.8838803 
##                coef. of skewness  =  -0.1231663 
##                coef. of kurtosis  =  2.542364 
## Filliben correlation coefficient  =  0.9962534 
## ******************************************************************

plot(gamlss(bron~poll_centered,data=Bronchitis,family=BI))

## GAMLSS-RS iteration 1: Global Deviance = 213.2586 
## GAMLSS-RS iteration 2: Global Deviance = 213.2586

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.01887182 
##                        variance   =  1.227365 
##                coef. of skewness  =  -0.1119096 
##                coef. of kurtosis  =  3.24346 
## Filliben correlation coefficient  =  0.9963054 
## ******************************************************************

plot(gamlss(bron~poll_centered,data=Bronchitis,family=BI))

## GAMLSS-RS iteration 1: Global Deviance = 213.2586 
## GAMLSS-RS iteration 2: Global Deviance = 213.2586

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.07470508 
##                        variance   =  1.137648 
##                coef. of skewness  =  -0.04981677 
##                coef. of kurtosis  =  2.853494 
## Filliben correlation coefficient  =  0.9982264 
## ******************************************************************

# plot fit
plot(Bronchitis$poll_centered,Bronchitis$bron,col="blue4",
        ylab = "Absence/Presense of Bronchitis", xlab = "Pollution level")
abline(h=c(0,1),col="light blue")
axe.x = seq(-10,10,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])/(1+exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2]))
lines(axe.x,f.x,col="pink2",lwd=2)

Model check suggests a good fit. Lets finally check if the interaction is significant:

# interaction ?
fit.glm = glm(bron~poll_centered*cigs,data=Bronchitis,family=binomial)
summary(fit.glm)

## 
## Call:
## glm(formula = bron ~ poll_centered * cigs, family = binomial, 
##     data = Bronchitis)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4389  -0.5698  -0.4144  -0.2927   2.3970  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -2.394978   0.294856  -8.123 4.57e-16 ***
## poll_centered       0.158671   0.071836   2.209   0.0272 *  
## cigs                0.214105   0.038699   5.533 3.16e-08 ***
## poll_centered:cigs -0.005292   0.010257  -0.516   0.6059    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 221.78  on 211  degrees of freedom
## Residual deviance: 173.95  on 208  degrees of freedom
## AIC: 181.95
## 
## Number of Fisher Scoring iterations: 5

anova(fit.glm,test="LRT")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: bron
## 
## Terms added sequentially (first to last)
## 
## 
##                    Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                                 211     221.78              
## poll_centered       1    8.519       210     213.26  0.003515 ** 
## cigs                1   39.044       209     174.21 4.143e-10 ***
## poll_centered:cigs  1    0.264       208     173.95  0.607347    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interaction is not significant.

(ii) myocardialinfarction.csv

The file myocardialinfarction.csv indicates if a participant had a myocardial infarction attack (variable infarction) as well the participant’s treatment (variable treatment).

Does Aspirin decrease the probability to have a myocardial infarction attack ?

Lets (i) import the dataset, (ii) change the levels of the factor treatment so that Placebo corresponds to the reference group, (iii) and finally plot the (sample) probabilities to get an attack by treatment group

# import
myocardialinfarction = read.csv("data/myocardialinfarction.csv")
# by default, Aspirin is the reference group as the alphabetic order is used
myocardialinfarction$treatment = factor(myocardialinfarction$treatment,
                                         levels=c("Placebo","Aspirin"))
# plot
par(mfrow=c(1,1),mar=c(3,4,3,1))
pi.group =  tapply(myocardialinfarction$infarction=="attack",myocardialinfarction$treatment,mean)
table.group = tapply(myocardialinfarction$infarction=="attack",myocardialinfarction$treatment,table)
temp = barplot(pi.group,plot=FALSE)
barplot(pi.group,ylab="Probability",xlab="",
        main = "Probability of myocardial infarction\n per treatment group",names=rep("",2),
        cex.axis=.6,axes=FALSE,cex.main=1.4)
axis(2,las=2,cex.axis=.8)
axis(1,temp[,1],names(pi.group),cex.axis=1.25,tick=FALSE)
for(gw in 1:length(pi.group)){
    text(temp[gw],pi.group[gw]/2,.p(table.group[[gw]]["TRUE"]," / ",sum(table.group[[gw]])),
         col="red",cex=1.5)
    }

The barplot seems to suggest that the treatment (aspirin) reduces the risk of myocardial infarction. Lets fit a logistic model to assess if this difference is significant. Note that, in this case (a dichotomous outcome and a dichotomous predictor), a test of equality of proportions or an independence test could also do the job. With a logistic model, other predictors could easily be added to the model and the beta parameter corresponding to the treatment can be interpreted by means of odd ratios (or relative risk ratios when prevalences are small, as we will note at the end of this practical).

# model fit
fit.glm = glm(I(infarction=="attack")~treatment,data=myocardialinfarction,family=binomial)
summary(fit.glm)

## 
## Call:
## glm(formula = I(infarction == "attack") ~ treatment, family = binomial, 
##     data = myocardialinfarction)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.1859  -0.1859  -0.1376  -0.1376   3.0544  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -4.04971    0.07337 -55.195  < 2e-16 ***
## treatmentAspirin -0.60544    0.12284  -4.929 8.28e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3114.7  on 22070  degrees of freedom
## Residual deviance: 3089.3  on 22069  degrees of freedom
## AIC: 3093.3
## 
## Number of Fisher Scoring iterations: 7

# test of equality of (independent) proportions
prop.test(unlist(lapply(table.group,function(x)x[2])),unlist(lapply(table.group,sum)))

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  unlist(lapply(table.group, function(x) x[2])) out of unlist(lapply(table.group, sum))
## X-squared = 24.429, df = 1, p-value = 7.71e-07
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.004597134 0.010814914
## sample estimates:
##     prop 1     prop 2 
## 0.01712887 0.00942285

# test of independence
chisq.test(matrix(unlist(table.group),ncol=2))

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  matrix(unlist(table.group), ncol = 2)
## X-squared = 24.429, df = 1, p-value = 7.71e-07

The three methods lead to the same conclusion: there is a significant difference between the probabilities of having a myocardial infarction of the two treatment groups. Note that the two last methods get exactly the same results (as they use the same test X-squared statsitic).

Lets define the fitted probabilities to get an attack and compare them to the sample probabilities: they should match:

# pi_placebo:
pi_placebo = exp(-4.04971)/(1+exp(-4.04971))
pi_placebo

## [1] 0.01712891

pi.group[1]

##    Placebo 
## 0.01712887

# pi_aspirin
pi_aspirin = exp( -4.04971-0.60544)/(1+exp( -4.04971-0.60544))
pi_aspirin

## [1] 0.009422852

pi.group[2]

##    Aspirin 
## 0.00942285

Finally, note that when prevalence are small, the exponential of the logistic regression corresponding to the treatment may also be interpreted as relative risk ratios. Indeed:

# interpreation of exp(beta1) when prevalences are small
pi_aspirin/pi_placebo

## [1] 0.5501137

exp(-0.60544)

## [1] 0.5458342

c(exp(-0.60544-qnorm(.975)*0.12284),exp(-0.60544+qnorm(.975)*0.12284))

## [1] 0.4290414 0.6944202

Thus, aspirin strongly reduces the risk of myocardial infarction.

(ii) crabs.csv

This data set is derived from Agresti (2007, Table 3.2, pp.76-77). It gives 6 variables for each of 173 female horseshoe crabs:

Explanatory variables that are thought to affect this included the female crab’s color (C), spine condition (S), weightweight (Wt)
C: the crab’s colour,
S: the crab’s spine condition,
Wt: the crab’s weight,
W: the crab’s carapace width,
Sa: the response outcome, i.e., the number of satellites.

Check if the width of female’s back can explain the number of satellites attached by fitting a Poisson regression model with width.

Lets import the datasset, fit a poisson loglinear model, plot the fit and perfom a model check :

crabs = read.csv("data/crab.csv",header=TRUE)

# plot:
plot(crabs$W,crabs$Sa,col="blue4",
        ylab = "Number of satellites", xlab = "width of female's back")
abline(h=c(0),col="light blue")

# fit
fit.glm = glm(Sa~W,data=crabs,family=poisson)
library(gamlss)


# plot fit:
plot(crabs$W,crabs$Sa,col="blue4",
        ylab = "Number of satellites", xlab = "width of female's back")
abline(h=c(0),col="light blue")

axe.x = seq(15,40,length=1000)
f.x = exp(fit.glm$coef[1]+axe.x*fit.glm$coef[2])
lines(axe.x,f.x,col="pink2",lwd=2)

    # not a great fit...

# model check
plot(gamlss(Sa~W,data=crabs,family=PO))

## GAMLSS-RS iteration 1: Global Deviance = 923.1762 
## GAMLSS-RS iteration 2: Global Deviance = 923.1762

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.1836972 
##                        variance   =  2.855978 
##                coef. of skewness  =  0.5993409 
##                coef. of kurtosis  =  2.564461 
## Filliben correlation coefficient  =  0.9717529 
## ******************************************************************

plot(gamlss(Sa~W,data=crabs,family=PO))

## GAMLSS-RS iteration 1: Global Deviance = 923.1762 
## GAMLSS-RS iteration 2: Global Deviance = 923.1762

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.2039895 
##                        variance   =  2.938738 
##                coef. of skewness  =  0.5165437 
##                coef. of kurtosis  =  2.603542 
## Filliben correlation coefficient  =  0.9819214 
## ******************************************************************

plot(gamlss(Sa~W,data=crabs,family=PO))

## GAMLSS-RS iteration 1: Global Deviance = 923.1762 
## GAMLSS-RS iteration 2: Global Deviance = 923.1762

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.1817619 
##                        variance   =  2.766722 
##                coef. of skewness  =  0.566886 
##                coef. of kurtosis  =  2.728251 
## Filliben correlation coefficient  =  0.9772948 
## ******************************************************************

    # confirm lack of fit -> bin the estimates
    
# 2 alternative models
plot(gamlss(Sa~W,data=crabs,family=ZIP))

## GAMLSS-RS iteration 1: Global Deviance = 766.424 
## GAMLSS-RS iteration 2: Global Deviance = 760.0696 
## GAMLSS-RS iteration 3: Global Deviance = 760.0688

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  0.03640043 
##                        variance   =  1.195172 
##                coef. of skewness  =  0.6257093 
##                coef. of kurtosis  =  4.205008 
## Filliben correlation coefficient  =  0.9849825 
## ******************************************************************

plot(gamlss(Sa~W,data=crabs,family=NBI))

## GAMLSS-RS iteration 1: Global Deviance = 751.2934 
## GAMLSS-RS iteration 2: Global Deviance = 751.291 
## GAMLSS-RS iteration 3: Global Deviance = 751.291

## ******************************************************************
##   Summary of the Randomised Quantile Residuals
##                            mean   =  -0.01039779 
##                        variance   =  1.043264 
##                coef. of skewness  =  -0.1133986 
##                coef. of kurtosis  =  2.291098 
## Filliben correlation coefficient  =  0.9929856 
## ******************************************************************

    # check ?ZIP and ?NBI for detail

Reasonably, there is a lack of fit> the estimates are not to be trusted.

GLM with R

D.-L. Couturier / R. Nicholls / M. Fernandes

Last modified: 21 Nov 2022

Section 1: Logistic regression

Section 1.1: importation and descriptive analysis

Section 1.2: Model fit

Section 1.3: Model selection

Section 1.3: Model check

Section 1.4: Fun

Section 2: Poisson regression

Section 2.2: Importation

Section 2.2: Model fit

Section 2.3: Model selection

Section 2.3: Model check

Section 6: Practicals

(i) Bronchitis.csv

(ii) myocardialinfarction.csv

(ii) crabs.csv